# Library
library('rjags')

#CHANGE WORKING DIRECTORY
setwd("/Users/felipenunes/Desktop/Pesquisas/15.Left presidentsRIEL/SUBMISSAO/Bayesian/")

# Data set
df <- read.table("data.csv", sep=",", header=T)
df$interac <- df$strong.party*df$X2round.hap
names(df)

# SEPARATE X & Y
y <- df[,2]
x <- as.matrix(df[,c(8, 10, 9, 21, 11, 23)])
head(x)

# DIFFERENT PRIORS TO TRY
# N is sample size
# K is number of parameters beta to estimate
# m is the prior mean for each parameter
# prec is the prior precision for each parameter
# mbeta0 is mean prior for beta0
# precbeta0 is precision prior for beta0
# y is depend. variable
# x is the matrix of independent variables
dataA <- list(n=length(y), K=6, m=c(2, 4, -1, -1, -1, 5), prec=c(.021, .011, .020, .013, .012, .034), x=x, y=y)

#SET UP INITAL VALUES
# This creates a list with 5 copies of the initial values.
# Change the number 5 to match the n.chains variable in the bugs call.
inits <- rep(list(list(beta=c(1,1,1,1,1,1))), 5)

# DEFINE PARAMETERS TO MONITOR
parameters <- c("beta")

## GLM Logit
m1 <- glm(left.type2 ~ X2round.hap*strong.party + oil.price + poverty + union.support + growth.before + vote.margin - 1, data=df, family=binomial(link=logit))
summary(m1)

## invoke JAGS
lab2.sim <- jags.model(file="logit.bug", data=dataA, n.chains = 5)

## Results
out <- coda.samples(lab2.sim, n.iter=2000, thin=10, variable.names="beta")

summary(out)

plot(out)

# Using the mean of all chains
temp1 <- out[[1]]
temp2 <- out[[2]]
temp3 <- out[[3]]
temp4 <- out[[4]]
temp5 <- out[[5]]

temp <-as.data.frame(rbind(temp1, temp2, temp3, temp4, temp5))
head(temp)

# confidence Intervals
mysummary = function(invector) {
  c(mean(invector), sd(invector), sd(invector)/sqrt(length(invector)), quantile(invector, .025), quantile(invector,.975))
}

tab <- data.frame(t(apply(temp2, 2, mysummary)))
colnames(tab) <- c('Coeff', 'SD', 'SE', '2.5%', '97.5%')
rownames(tab) <- c('Run-off', 'Strong Party', 'Oil price', 'Poverty', 'Union support', 'Interaction')
tab

# Parameters estimated with MCMC
round2 <- temp[,1]
party <- temp[,2]
oil <- temp[,3]
pov <- temp[,4]
union <- temp[,5]
inter <- temp2[,6]

# Draw Kernel Density Plots for each parameter distribution. The densityplot methods and qqmath methods produce empirical density and probability plots.
densityplot(out)

# Plotting Auto-correlation of estimates in the series 
par(mfrow=c(3,2))
acf(round2, plot=T)
acf(party, plot=T)
acf(oil, plot=T)
acf(pov, plot=T)
acf(union, plot=T)
acf(inter, plot=T)

# Traceplot
par(mfrow=c(3,2))
traceplot(out)

# Scatterplot matrix of correlation plots
splom(temp[ ,1:6])

# Other checkings
geweke.diag(out)
heidel.diag(out)
raftery.diag(out)
effectiveSize(out)

## Canned Bayesian Package
require(MCMCglmm)
sigma2 <- 1
k <- sqrt(1 + (((16 * sqrt(3))/(15 * pi))^2) * sigma2)

m <- MCMCglmm(left.type2 ~ X2round.hap + strong.party + oil.price + poverty + union.support + growth.before + vote.margin - 1,
  data=na.omit(df), family = "categorical", prior = list(
    B = list(mu = rep(0, 7), V = diag(7) * 100),
    R = list(V = sigma2, fix = 1)),
  nitt = 1000 * 5000 + 20000, thin = 5000, burnin = 20000, verbose=FALSE)

autocorr(m$Sol)
plot(m$Sol)

## regular model
summary(m)

## MCMCglmm is on a different scale than usual logit, rescale by k
## to put on regular logistic scale
cbind(M = colMeans(m$Sol), HPDinterval(m$Sol))/k